import os

os.system("echo q | gmx make_ndx -f md_0_1.gro")

###分析
#周期性边界
os.system("echo 1 | gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC_MOL.xtc -pbc nojump -n index.ndx")

# MMGBSA
#os.system("cp -rf ../amber99sb-ildn.ff ./")
#os.system("conda activate gmxMMPBSA && screen -L -dmS MMPBSA gmx_MMPBSA -O -i mmgbsa-LIG.in -ci index.ndx -cs md_0_1.tpr -cp topol.top -ct md_0_1_noPBC_MOL.xtc -cg 1 13")
#conda activate gmxMMPBSA
#screen -L -dmS MMPBSA gmx_MMPBSA -O -i mmgbsa-LIG.in -ci index.ndx -cs md_0_1.tpr -cp topol.top -ct md_0_1_noPBC_MOL.xtc -cg 1 13

#rmsd
os.system("echo 1 1| gmx rms -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o rmsd-pro.xvg -tu ns")
#os.system("echo 13 13| gmx rms -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o rmsd-lig.xvg")
#os.system("echo 20 20| gmx rms -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -n index.ndx -o rmsd-complex.xvg")
#回旋半径
os.system("echo 1 | gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o gyrate-pro.xvg")
#rmsf
os.system("echo 1 | gmx rmsf -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o rmsf-pro.xvg -oq bfactor.pdb -res")
#os.system("echo 20 | gmx rmsf -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc  -o rmsf-complex.xvg -oq bfactor-complex.pdb -res  -n index.ndx")
#H bonds
os.system("echo 1 1 | gmx hbond -f md_0_1.xtc -s md_0_1.tpr -n index.ndx -num hbond-pro-lig.xvg -hbn -hbm -tu ns")
#抽帧
os.system("echo 1 | gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o traj.pdb -dt 1000 -n index.ndx")
#SASA
os.system("echo 1 | gmx sasa -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o sasa-pro.xvg -tu ns")


#python LandScape-RMSD-Gyrate.py 
# print("protein-protein")
os.system("echo 1  1 | gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC_MOL.xtc -o md_0_1_noPBC_fit.xtc -fit rot+trans")
# 选择backbone进行计算和输出
os.system("echo 3  3| gmx rms -s md_0_1.tpr -f md_0_1_noPBC_fit.xtc -o rmsd-landscape.xvg")
# 同样选择backbone
os.system("echo 3 | gmx gyrate -s md_0_1.tpr -f md_0_1_noPBC_fit.xtc -o gyrate-landscape.xvg")



rmsd = open("rmsd-landscape.xvg","r")
rmsd_lines = [i for i in rmsd]
rmsd.close()

gyrate = open("gyrate-landscape.xvg","r")
gyrate_lines = [i for i in gyrate]
gyrate.close()


head = []
pc1_list = []
pc2_list = []
for i in range(len(rmsd_lines)):
    if rmsd_lines[i][0] == "@":
        head.append(rmsd_lines[i])
    else:
        pc1_list.append(rmsd_lines[i][0:-2])
# print(pc1_list)
for i in range(len(gyrate_lines)):
    if gyrate_lines[i][0] != "@":
        pc2_list.append(gyrate_lines[i][15:22])
# print(pc2_list)
for i in range(len(pc1_list)):
    pc1_list[i] = pc1_list[i] + "    " + pc2_list[i] + "\n"

# print(pc1_list)

pc12 = open("rmsd_gyrate_sham.xvg","w+")
pc12.writelines(head+pc1_list)
pc12.close()


os.system("gmx sham -tsham 310 -nlevels 100 -f rmsd_gyrate_sham.xvg -ls rmsd_gyrate_gibbs.xpm -g rmsd_gyrate.log -lsh enthalpy.xpm -lss entropy.xpm")
#os.system("gmx sham -tsham 300 -nlevels 100 -f pc12_sham.xvg -ls pc12_gibbs.xpm -g pc_12.log -lsh pc12_enthalpy.xpm -lss pc12_entropy.xpm")

# os.system("pip install DuIvyTools  -i https://pypi.tuna.tsinghua.edu.cn/simple")
# os.system("dit xpm_show -f pc12_gibbs.xpm -ip -x RMSD -y Gyrate")
# os.system("dit xpm_show -f pc12_gibbs.xpm -3d -x RMSD -y Gyrate")
os.system("dit xpm_show -f rmsd_gyrate_gibbs.xpm -ip -x RMSD -y Gyrate -o LandScape-2d.tif -ns")
os.system("dit xpm_show -f rmsd_gyrate_gibbs.xpm -3d -x RMSD -y Gyrate -o LandScape-3d.tif -ns")

#collaction
os.system("mkdir download")
os.system("cp LandScape-*tif ./download")
os.system("cp *-pro.xvg ./download")
os.system("cp *-lig.xvg ./download")
os.system("cp *-complex.xvg ./download")
os.system("cp bfactor* ./download")
os.system("cp traj.pdb ./download")
os.system("cp bindex.ndx ./download")
os.system("cp rmsd_gyrate_gibbs.xpm ./download")
os.system("cp rmsd_gyrate.log ./download")

